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ABSTRACT 

We investigate hydrodynamic instability of accelerating ionization fronts us- 
ing two dimensional hydrodynamic simulations that include detailed energy de- 
position and release due to the absorption of UV radiation, recombination of 
hydrogen, radiative molecular cooling, and magnetic pressure. We consider lin- 
ear perturbation growth and find that the stabilization mechanism associated 
with non-accelerated fronts remains a significant factor even when acceleration is 
present. Conversely, if recombination in the ionized region is turned off, Rayleigh- 
Taylor (RT) instability becomes effective, and the classical RT growth rate re- 
covered. 

Subject headings: HII regions - ISM: molecules - ISM: kinematics and dynamics 
- hydrodynamics - instabilities - methods: numerical 



1. INTRODUCTION 

Ionization fronts (IF) near OB stars have attracted considerable interest due to their 
amazing shapes. For example, the Eagle Nebula has three famous 'pillars' or 'elephant 
trunks' (Hester et al. 1996; Pound 1998). Strong UV radiation from nearby OB stars illumi- 
nates the surface of a molecular cloud, and photoevaporation occurs in a very thin layer at 
the surface. As a result, an ablation flow begins and a shock propagates into the cloud. Such 
an ionization or ablation front is categorized as D-type. The size of the system including the 
OB stars and the molecular cloud is a few parsecs, and the width of pillars is usually less 
than a parsec. 
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A number of previous theoretical works and numerical simulations (Williams, Ward- 
Thompson, & Whitworth (2001), Williams (2002) and references therein) have addressed 
understand the formation of pillars, but this problem is still an area of active research. All 
the models considered thus far can be divided into two large categories: those that relate 
formation of pillars to the presence of pre-existing dense clumps 'excavated' by moving 
ablation front, and those that relate the pillars to a nonlinear stages of ionization front 
instabilities. In this article we study the second category. The first category has been 
also extensively studied by various authors, most recently by Williams, Ward- Thompson, & 
Whitworth (2001), who reference earlier work. 

In modeling ionization front instability, researchers have pursued one of two approaches: 
they have either studied instability of the RT type, associated with finite acceleration of 
the cloud driven by the ablation pressure, or studied instability of a non-accelerating front 
associated with details of ionization/recombination processes in the ionized outflow. As far as 
we know, our paper will be the first where both the acceleration and ionization/recombination 
processes are taken into account simultaneously. 

To put our paper into historical perspective, we mention that RT instability of molecular 
clouds was studied by Spitzer (1954) and Prieman (1954). The case where the incident 
radiation is tilted with respect to the cloud surface was studied in Ryutov ct al. (2003). 
Instability of non- accelerated ionization fronts was studied by Vandervoort (1962), with the 
radiation tilt included. Kahn (1958) argued qualitatively that absorption of the incoming 
radiation by the hydrogen atoms formed by recombination in the ablation outflow can have 
a stabilizing effect, because of the stronger absorption near the dimples of the surface relief 
compared to the bumps. Axford (1964) presented a quantitative study that showed that this 
stabilization mechanism is most effective for perturbations with wavelengths larger than the 
recombination length. Sysoev (1997) provided more complete analysis and find the growth 
of long wavelength instabilities for normally incident radiation. Williams (2002) confirmed 
it and included effects of the radiation tilt. 

The main result of our study is that the linear RT instability, which is usually very 
robust and hard to stabilize, is actually also stabilized by the recombination effect. 

2. BASIC EQUATIONS 

Most theoretical analyses and numerical simulations assume that the gas is isothermal 
and do not take into account the energy budget. However it is possible to formulate the 
problem including detailed energy treatment. At the ablation front, an incident photon 
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whose energy is higher than the Lyman hmit is absorbed, ionizing the neutral gas and 
increasing its thermal energy. In the ionized region ionized and electrons recombine and 
thermal energy is radiated . When an ablatively driven shock propagates into the molecular 
gas, the gas is radiatively cooled on a very short time scale (Ncufcld, Lepp, & Melnick 
1995). The temperature of the interior of the molecular cloud in the Eagle Nebula pillars is 
a few tens of Kelvins (Pound 1998). In our work we solve the equations of two dimensional 
hydrodynamics with energy source terms that account for the radiative and recombinative 
processes just described: 

|^ + V-(pn) = 0, (1) 
^ + V-(pwtz + pI) = 0, (2) 

^ ('^ (^"^ ^ ^) ) ^ ^ ( ('^ (^^^ + + = -qre + quv - Qmol- (3) 

In order, these are the mass, momentum and energy conservation equations, p is mass 
density, p is pressure, I is unit tensor, u is the velocity vector, and e is the specific internal 
energy, g^e, Quv, and q„ioi are energy source terms due to the recombination in the ionized 
region, absorption of UV radiation from OB stars, and cooling in the molecular gas due to 
rotation and vibration mode of hydorgen molecules, respectively. 

Photoionization and recombination are considered with the following equations: 

df 

n— + nu-\^f = an{l-f)J-aBn^f, (4) 

^^-an{l-f)J, (5) 

where / = rii/n is the ionization fraction (/ = or 1 correspond to neutral or fully ion- 
ized gas, respectively), n is the total number of hydrogen nuclei per unit volume, including 
nuclei in molecules, atoms and ions, and rii is ionized hydrogen volume density. For numer- 
ically simplify, we do not take into account for the state of the neutral hydrogen, assuming 
that all molecular hydrogen dissociates before the IF, although the photodissociation region 
dominates some thickness (Tielens and HoUenbach 1985) which is longer than that of pho- 
toionization front. The temperature in the HI region becomes about a few tens to a few 
hundreds Kelvin. Our treatment is just as same as the isothermal models done in a number 
of theoretical and numerical works so far. We also do not consider the energy budget of the 
dissociation process, a — 6x 10~^^cm^ is photoionization cross-section of hydrogen, and J is 
the number flux of ionizing photons whose energy exceeds the critical energy for the ioniza- 
tion of hydrogen (13.6eV). We assume that the incident photons are parallel to the y axis. 
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We do not include recombination to the ground state, instead assuming diffuse radiation and 
its absorption are balanced locally (on the spot approximation), since the effect of this diffuse 
photon radiation is estimated at 10% ~ 20% (Canto et al. 1998). Only the so-called "case B" 
recombination {as = 2.6 x 10~^^cm^s~^ at T = lO'^K (Hummer & Scaton 1963)) is considered 
and fixed in this study, where as is the sum of the recombination coefficients to all hydro- 
gen states except to the ground state. Energy sources are written as: g^e = (nf'Y PskBT , 
q^^ = Wna{l - f)J, q^oi = n^^^ x IQ-^^erg cm-^s-^, and T = {rrip/h) x [4e/(7/ + 5)], 
rijnoi = ^^-(1 — /)/2, where T is temperature in Kelvin (pe = 2 x 1.5{nf)kBT + 2.5njnoikBT), 
nip is the proton mass, kj, is the Boltzmann constant, and rimod is molecular hydrogen vol- 
ume density. We use {nfyPskbT instead of {3/2){nfyaBkf,T for the recombination cooling 
term to model the thermal velocity dependence of the rates of recombination and free-free 
collisional cooling. /Sb = 1.25«b at T = 10^ K (Hummer & Seaton 1963) is also fixed. It 
is assumed that the averaged energy of the incident photon is (13.6 + W) eV per photon 
from O stars. The energy W is deposited into the gas as internal energy when a neutral 
hydrogen atom absorbs an incident photon. The ionized region becomes isothermal quickly 
due to the energy balance between cooling and heating processes. We take W — 1.73 x 10~^^ 
erg to produce an isothermal temperature of T = 10^ K in ionized gas. The cooling term 
Qmoi is effective for 40K <T < 3000K. Since this coohng effect is enough strong to keep the 
temperature in the molecular gas 40 K, this treatment is very similar to isothermal model. 
The more realistic treatment was done by Richling & Yorke (2000) for the calculations of 
the photoevaporation of protostellar disks. The main difference is the introduction of a 
magnetic pressure described soon. We ignore the metal line cooling for simplicity, since the 
cooling power has same dependence on the recombination cooling power (oc {nfY) assuming 
an equilibrium temperature. Our conclusion mainly does not change, whether we take into 
account the metal line cooling or not. The temperature T in the ionized outflow can be 
considered as a free parameter as it is controlled by the parameter W . In this respect, an 
explicit accounting for the radiative cooling does not seem to be critical to our analysis. We 
should note that the heating function W depends on the temperature of the surface of OB 
stars. 

To close the equations we also solve the equation of state: 



The first term on the right hand side is thermal pressure (/ = or 1 give (2/5)pe or (2/3)pe 
corresponding to diatomic or monoatomic adiabatic gas) and the second one is magnetic 
pressure, p^^ and are constant values. The index 7m (= 4/3 in this study) is also 
constant. can be 4/3 for magnetic turbulence or 2 for an initially uniform magnetic field. 
This magnetic pressure is introduced to prevent radiative collapse due to molecular cooling 
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when the shock heating occurs (Ryutov & Remington 2002). When the gas is compressed, 
the second term acts as a stiffness. 

The sound speed is then 

The hydrodynamic equations are solved with a Godunov-type code. The numerical 
scheme is the same as used in Mizuta, Yamada, & Takabe (2002), extended to accept real 
gas equation of state (Glaister 1999). Photoionization and recombination in Eq. 4 and 5 are 
solved implicitly at both half and full step (WiUiams 1999). 

3. NUMERICAL CONDITIONS 

We use two dimensional plane geometry (x-y). The grid size is uniform {Ax = Ay = 
2.5 X lO^^pc ~ 7.5 X lO^^cm which is still longer than the mean free path of an incident 
photon for the molecular cloud). 

Initially a 0.25 pc thick cloud is set 0.5pc away from the boundary {y = 3 pc) where 
the incident photons are introduced. The cloud has initial hydrogen number density n(H) = 
lO^cm"^ and temperature T = 40 K. Ionized gas is put between the boundary and the cloud 
surface. The number density is n(H)= lOcm"'^ and the gas is pressure matched with the 
molecular cloud. Behind the cloud a dilute molecular gas is put (n(H)= 10cm~^ and T = 40 
K) . Then this rear cloud surface is not pressure matched initially. Constant parameters 
and introduced in our EOS are chosen to be the thermal pressure and mass density, 
respectively, of the initial molecular cloud. The photon flux is parallel to the y-axis and 
constant except for a short imprinting period. All of the gases are at rest at t — when the 
radiative flux is turned on. 

The dynamics is divided into three phases. At first, an ablation flow begins and a shock 
propagates into the molecular cloud (I). Since the gas is assumed to be at rest initially, the 
ablation front moves in this frame. This is more realistic than simulating in the frame where 
the ablation front is at rest because we can include the effect of increasing distance of the 
cloud surface from the radiation source. This effect becomes important in the case with 
recombination. The velocity of the ablation front is almost constant in phase I. When the 
shock breaks out from the back side of the cloud, a rarefaction appears and proceeds to the 
ablation front (II). This phase is very short, less than lOky in most cases, although its length 
depends on incident photon intensity, number density of the cloud, and cloud thickness. 
When the rarefaction arrives at the ablation front, the acceleration phase begins (III). 
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A perturbation is introduced into the incident photon number flux for 10 ky when phase 
III begins. The perturbation is sinusoidal with a 10% amphtude J = Jo(l + 0.1cos(27ra;/A)), 
where Jq is the unperturbed incident photon number flux and A is the wave length. We 
studied three cases with different wave lengths: A = 0.92, 0.6, and 0.46 pc, which are similar 
to the width of the pillars of the Eagle nebula. We employ periodic boundary conditions 
at both edges that are parallel to y the axis, and free boundary conditions at the other 
boundaries. We follow the evolution at least until the ablation front propagates a distance 
equal to the wavelength of the perturbation. This means the distance between the cloud 
surface and the boundary where the photon flux is incident becomes longer than 1.5 pc in 
later times. 



4. RESULTS AND DISCUSSION 

4.1. Without Recombination 

At first, we discuss the case without recombination, namely, as — Pb — ^- The incident 
photon fiux at the boundary is | Jo| = 2.6 x 10^cm~^s~^. Figure 1 shows the incident photon 
number flux (solid lines) and the (/ = 0.5) ablation front contour at t = 340 ky. The 
absorption occurs in a very thin layer at the ablation front because of the very short photon 
mean free path because there is no absorption in the ionized gas. Figure 2 shows the time 
evolution of the perturbation amplitude, where the amplitude is deflned as the half distance 
between the maximum and minimum points in y of the contour / = 0.5. Results are 
shown for the three wavelengths. The perturbation grows with time for all wavelengths, and 
the growth rate is in good agreement with classical RT theory (Chandrasekhar 1981) (the 
amplitude is A = Aoexp(7t), where 7 = {kgY^"^ is the growth rate, k — 27r/A is the wave 
number, and g is the effective gravity, determined from a one dimensional simulation without 
any perturbation. The atwood number is assumed to be unity.) In Figure 2, three short 
segments shown shown used in the comparison to RT theory. Without recombination, the 
very small imprinted perturbation (comparable to the grid size), produces non-linear growth 
at later times. 



4.2. With Recombination 

On the contrary, in the case with recombination, the behavior of the imprinted perturba- 
tion is quite different. The incident photon number flux is increased to | Jo| = 5 x 10^^cm~^s~^ 
so that the effective gravity in the acceleration phase is almost the same as in the case without 



- 7- 



recombination. This effective flux is possible for the Eagle Nebula, where the total number 
flux from O stars is 2 x 10^°s~^ (Hester et al. 1996) and the temperature near the IF is 10^ K, 
although it depends on the distance from the radiation source, and on the density structure 
and temperature of the ionized gas. With this intensity and initial molecular cloud density, 
the number density of the ionized gas at the ablation front becomes about lO^cm"^. The 
hydrogen number density of shocked molecular gas is ~ 3 x lO^cm"^. These densities are 
consistent with observations (Hester et al. 1996; Pound 1998). Figure 2 shows time evolu- 
tion of the amplitude for the three wavelengths, with recombination. The number density 
of ionized gas at the ablation front decreases when the photon number flux decreases. The 
perturbation does not grow, contrary to the case without recombination, and the amplitude 
is oscillating in time. 

The qualitative explanation of the absence of the growth in this case is very similar to 
that offered by Kahn (1958) and Axford (1964). Some of the ablated material in the bubble 
(concave) region converges and a dense region appears in the ionized gas near the ablation 
front (Fig. 3). This density structure affects the absorption near the ablation front (Fig. 
4). In the dense gas the recombination rate, which has a quadratic dependence on density 
(Eq. 4) is high, and the absorption of incident photons is enhanced. The photon flux is 
larger around the spike, which is strongly driven so that the rocket effect is strong there. 
On the other hand, incident photons are absorbed only moderately around the bubble, and 
the rocket effect is weaker there. As a result, the amphtude is reduced. Once the amphtude 
is reduced multi-mode patterns appear and oscillate but do not grow, contrary to the case 
without recombination. 

This large difference between the cases with and without recombination results from 
a competition between RT push (wanting to make perturbations grow), and recombina- 
tion/opacity damping (tending to make perturbations anneal and flatten out). Only the RT 
push works without recombination. On the other hand with recombination both effects are 
available. Then, a small amplitude perturbation does not grow because the damping works 
effectively. This suggests there may be a critical amplitude for the RT push to win and pro- 
duce nonlinear growth. The growth rate is only weakly affected by ablative stabilization (for 
example Takabe, Mima, Montierth, & Morse (1985)), because the wavelengths of interest 
are quite large; as we consider only normal incidence of the radiation, effects of radiation tilt 
Ryutov et al. (2003) are also absent. 

We also studied the same problem using Williams's isothermal model Williams (1999) for 
comparison. He assumed that both the ionized and molecular gases are isothermal instead of 
including a detailed energy budget accounting for radiative processes. The isothermal sound 
velocity is defined as = (1 -|- 99/)^/^ x 1 km s~^. The pressure and specific internal energy 
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are reset using this sound velocity and an adiabatic equation of state. The adiabatic index 
is assumed to be F = 1.1. The results are very similar to what we have shown here both 
without and with recombination, namely, the perturbation grows without recombination and 
does not grow and oscillates with recombination. Hence, this stabilization effect in the linear 
regime is observed using two different methods of simulating the radiatively driven cloud. 

5. CONCLUSION 

Instabilities at ionization fronts are discussed with effective gravitation using parameters 
seen in the Eagle Nebula. The acceleration phase has not been studied well and may occur 
in the evolution of Nebula. The perturbation on the ionization front grows in the case 
without recombination. The growth rate is in good agreement with classical RT instability. 
When recombination is turned on, which is the more realistic case, the difference of density 
profile causes a different absorption profile. This works to effectively smooth the surface. 
The perturbation does not grow and the amplitude oscillates in time. This damping effect 
works because with recombination the ablated plasma is not optically thin. We observed 
the same results using Williams' (1999) isothermal model. Linear regime perturbations 
appear to exhibit a stabilization effect in these IF simulations. In this study, the imprinted 
perturbation is very small and we observed only the liner growth regime with recombination. 
The dynamics may be quite different in the nonlinear regime. There may be a critical angle 
of the surface perturbation above which the RT push wins out in the competition between 
RT push and recombination/opacity damping. 
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Fig. 1. — Photon number flux contour (solid lines), and ablation front (/ = 0.5) contour 
(dashed line) at t = 340ky (without recombination case). Incident photon flux is incoming 
from the top of the flgure (at y=3 pc). Photon number flux contours are drawn every 
1.3 X 10®cm~^s~^ in linear scale at equal intervals. However, most contours override each 
other because of the very short mean free path of the incident photon in the molecular cloud. 
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Fig. 2. — Time evolution of log scale amplitude of perturbation in each wavelength (without 
recombination). Three wavelengths are shown: A = 0.92 pc (solid), A = 0.6 pc (dotted), 
A = 0.46 pc (dot dash). The short wavelengths grow faster than the long wavelengths, 
in accordance with classical RT theory. The short line segments shown are used for the 
comparison to classical RT growth. Results with recombination cases are also shown: A = 
0.92 pc (long dased), A = 0.6 pc (short dashed), A = 0.46 pc (two-dot dash). The shorter 
wavelengths appear to produce damped oscillation, whereas the longest wavelength produces 
oscillation. 
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Fig. 3. — Log scale of hydrogen number density contours at t = 150 ky for the simulation 
including recombination. Contours are drawn every 0.05 in log scale at equal intervals. 
Number density decreases away from the ablation front (see Fig.5). 
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Fig. 4. — Photon number flux and ablation front contours (similar to Fig.l) for simulations 
including recombination at the same time as the results shown in Fig.4. Contours are drawn 
every 2.5 x 10^° cm~^s~^ in linear scale at equal intervals. Flux decreases from top to bottom. 



